% /////////////////////////////////////////////////////////////////////////
% Paper: The Technical Default Spread
% Authors: Emilio Bisetti, Kai Li, and Jun Yu
% Date: 2024, May 14
% This code is used to conduct firm level simulation
%% //////////////////////////////////////////////////////////////////////////

clear; close;clc; 
% choose the model with risk shocks
cd([cd '/Model with Risk Shocks']);

% choose the model with covenant
cd([cd '/GE Model computation with covenant']);

load para_v1.mat;
load('mainmodel_tfpunc_results.mat'); %% save dynare simulation results

% You need to run the dynare code to get following result
% dynare mainmodel_tfpunc.mod     


% freq = 4;
num_endo = length(M_.endo_names); % Set the length of endogenous variables
endo_list = cell(num_endo, 1); % Set a vector to store the string of variables' names
simulated_data = oo_.endo_simul(:, end-1200+1:end);

for idx = 1:length(M_.endo_names)
    endo_list{idx} = M_.endo_names(idx,:);
end

for idx = 1:length(endo_list)
%     disp(i)
    eval([endo_list{idx}, '= transpose(simulated_data(idx,:));']);
end

% idiosyncratic shock Z_it
T = length(dy); % length of time period
N = 5000;       % number of firms
kappa = 1-eta; % the cutoff: entrpreneur below it get liquidation shock


% specify a starting value for F (in level) 
% recover the K with trend
length_Q = length(dk); % Length of quarters
K_Q = zeros(length_Q,1); 
K_Q(1) = 1; % assign initial value to be 1
for index = 2:length(dk)
   K_Q(index) = K_Q(index-1)*exp(dk(index));  % recover the time series of K_t
end

% These are normalized variables
Yt = exp(y).*K_Q;
Ct = exp(c).*K_Q;
It = exp(ii).*K_Q;
Nt = exp(n).*K_Q;
Bt = exp(b).*K_Q;
dY_Q = diff(Yt)./Yt(1:end-1);

% For prices that are not normalized
Qt = exp(q);
Ht = exp(h);
RKt = exp(rk);
RDt = exp(rd);
RLt = exp(rl);
act_deft = act_def; % actual default probability in level
TDprobt = TDprob; % technical default probability
sigma0t = exp(sigma0); % Time varing sigma0
sigma1t = exp(sigma0); % Time varing sigma0
omega0bart = omega0bar;

% initialize a panel of random shocks, 
% generate three set of random variables: omega0, omega1, and survival index

rng('default');
%omega1   = normrnd(mu_1, sig1, [T,N]);
surv_ind = rand(T,N); % Uniformly distributed random numbers

count = 0;

% Firm simulation
% one state variables 
Ni = nan*ones(T,N);
Ni(1,:) = Nt(1); % use the same value as initial value

% return
ret_panel = nan*ones(T,N);
% Stictness
Strict_panel = nan*ones(T,N);
Strict_panel2 = nan*ones(T,N);
K_panel = nan*ones(T,N);
omega1_hat= nan*ones(T,N);
omega0 = nan*ones(T,N);
omega1= nan*ones(T,N);

Signal_panel = nan*ones(T,N);
noise_panel = nan*ones(T,N);
noise_para = 1;
noise_sig=0.2;
TD_indicator = nan*ones(T,N);
max_signal=nan*ones(T,1);

for t_idx = 1:T-1  % time index 
    
    rng(t_idx); % set seed for each period
    omega0(t_idx+1,:) = normrnd(-0.5*sigma0t(t_idx)^2, sigma0t(t_idx), [1,N]);
    omega1(t_idx+1,:) = normrnd(-0.5*sigma1t(t_idx)^2, sigma1t(t_idx), [1,N]);
    
    
    for n_idx = 1:N  % firm index
        % at time t_idx +1, use H and RL, omegabar at t_idx, but use RK at
        % time t_idx+1
         if omega0(t_idx+1,n_idx) >= omega0bart(t_idx)
            noise_panel(t_idx+1,n_idx) =normrnd(0, noise_sig,1);
            TD_indicator(t_idx+1,n_idx) = 0;
            
         elseif omega0(t_idx+1,n_idx) < omega0bart(t_idx)%+0.05*abs(omega0bart(t_idx))
             noise_panel(t_idx+1,n_idx)=-10000000;
             TD_indicator(t_idx+1,n_idx)=1;
             
         end
        
        Signal_panel(t_idx+1,n_idx) = omega0(t_idx+1,n_idx)+noise_panel(t_idx+1,n_idx);
        Strict_panel(t_idx+1,n_idx) = -(omega0(t_idx+1,n_idx)-omega0bart(t_idx));
        omega1_hat(t_idx+1,n_idx) = h(t_idx)+rl(t_idx)-rk(t_idx+1)-log(1+exp(h(t_idx)))-omega0(t_idx+1,n_idx);
        
        if surv_ind(t_idx+1,n_idx) > kappa % those who are not hit by liquidation shock
              
            % get firm level net worth  
            if omega0(t_idx+1,n_idx) >= omega0bart(t_idx) && omega1(t_idx+1,n_idx) >= omega1_hat(t_idx+1,n_idx) % entrepreneur in control and no default
                Ni(t_idx+1,n_idx) = (exp(omega0(t_idx+1,n_idx)+omega1(t_idx+1,n_idx))*exp(rk(t_idx+1))*(1+exp(h(t_idx)))-exp(rl(t_idx))*exp(h(t_idx)))*Ni(t_idx,n_idx);
                ret_panel(t_idx+1, n_idx) = Ni(t_idx+1,n_idx)/Ni(t_idx, n_idx);
                
            elseif omega0(t_idx+1,n_idx) >= omega0bar(t_idx) && omega1(t_idx+1,n_idx) < omega1_hat(t_idx+1,n_idx)
                Ni(t_idx+1,n_idx) = chi*Nt(t_idx);
                ret_panel(t_idx+1, n_idx) = 0;                                
            else 

                Ni(t_idx+1,n_idx) = exp(rd(t_idx))*Ni(t_idx, n_idx);
                ret_panel(t_idx+1, n_idx) = Ni(t_idx+1,n_idx)/Ni(t_idx, n_idx);
                
            end
                

        elseif surv_ind(t_idx+1,n_idx) <= kappa
             
            % get firm level net worth  
             if omega0(t_idx+1,n_idx) >= omega0bart(t_idx) && omega1(t_idx+1,n_idx) >= omega1_hat(t_idx+1,n_idx) % entrepreneur in control and no default

                Ni(t_idx+1,n_idx) = (exp(omega0(t_idx+1,n_idx)+omega1(t_idx+1,n_idx))*exp(rk(t_idx+1))*(1+exp(h(t_idx)))-exp(rl(t_idx))*exp(h(t_idx)))*Ni(t_idx,n_idx);
                ret_panel(t_idx+1, n_idx) = Ni(t_idx+1,n_idx)/Ni(t_idx, n_idx);
                
            elseif omega0(t_idx+1,n_idx) >= omega0bar(t_idx) && omega1(t_idx+1,n_idx) < omega1_hat(t_idx+1,n_idx)

                Ni(t_idx+1,n_idx) = chi*Nt(t_idx);
                ret_panel(t_idx+1, n_idx) = 0;
                                
             else 
                Ni(t_idx+1,n_idx) = exp(rd(t_idx))*Ni(t_idx, n_idx);
                ret_panel(t_idx+1, n_idx) = Ni(t_idx+1,n_idx)/Ni(t_idx, n_idx);
                
            end
            
                 Ni(t_idx+1,n_idx) = chi*Nt(t_idx);
                 %ret_panel(t_idx+1, n_idx)= Ni(t_idx+1,n_idx)/Ni(t_idx, n_idx);
        end
        
    end
    
    % compute maximum value 
    non_td_indicator = omega0(t_idx+1,:)>= omega0bart(t_idx);
    non_td_omega0 = omega0(t_idx+1,non_td_indicator);
    non_td_signal = Signal_panel(t_idx+1,non_td_indicator);
    
    max_signal(t_idx+1) = max(omega0bart(t_idx)-non_td_signal);
    
    for n_idx = 1:N  % firm index
        % at time t_idx +1, use H and RL, omegabar at t_idx, but use RK at
        % time t_idx+1
         if omega0(t_idx+1,n_idx) >= omega0bart(t_idx)
            Strict_panel2(t_idx+1,n_idx) = 100*exp(omega0bart(t_idx)-Signal_panel(t_idx+1,n_idx))/(exp(max_signal(t_idx+1)));
         elseif omega0(t_idx+1,n_idx) < omega0bart(t_idx)%+0.05*abs(omega0bart(t_idx))
             Strict_panel2(t_idx+1,n_idx) = 100;
         end
    end
    
    
end

% We need to winsorize the return in the cross section each period
ret_panel_alt = ret_panel(2:end,:);
RET = nan(size(ret_panel_alt));
T_Ret = size(ret_panel_alt,1);
for j =1:T_Ret
    [RET_win,varargout] = winsor(ret_panel_alt(j,:)',[1 99]);
    RET(j,:)=RET_win';
end

size_panel = Ni(1:end-1,:);  
Strictness = Strict_panel2(2:end,:);
Signal = Strict_panel2(2:end,:);
TD_indicator = TD_indicator(2:end,:);

% sort portfolios�� compute size weighted portfolio level returns
[phy5] = portSort5VWQtr(Strictness', size_panel', RET');
Tuse = size(phy5,2);
exc_ret_port = phy5'-RDt(1:Tuse);
Ret_result = mean(exc_ret_port)*4*100;

% sort portfolios�� compute size weighted portfolio level returns
equal_weight = ones(size(size_panel'));
[strict5] = portSort5VWQtr(Strictness', equal_weight, Strictness');
aa = strict5';
strict5_result = mean(aa);

% report results of Panel B of Table 3
PanelB = [strict5_result;Ret_result]

